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We have performed a Quantum Monte Carlo study of a two-dimensional bulk sample of inter- 
acting 1/2-spin structureless fermions, a model of '^He adsorbed on a variety of preplated graphite 
substrates. We have computed the equation of state and the polarization energy using both the 
standard fixed-node approximate technique and a formally exact methodology, relying on bosonic 
imaginary-time correlation functions of operators suitably chosen in order to extract fermionic en- 
ergies. As the density increases, the fixed-node approximation predicts a transition to an itinerant 
ferromagnetic fluid, whereas the unbiased methodology indicates that the paramagnetic fluid is the 
stable phase until crystallization takes place. We find that two-dimensional '^He at zero temperature 
crystallizes from the paramagnetic fluid at a density of 0.061 with a narrow coexistence region 
of about 0.002 A~^. Remarkably, the spin susceptibility turns out in very good agreement with 
experiments. 



PACS numbers: 



I. INTRODUCTION 



A quasi-two-dimensional (2d) bulk ^He sample at zero 
temperature is a very fascinating scenario to explore the 
physics of strongly correlated Fermions. The liquid phase 
can be experimentally realized over a wide range of den- 
sities by adsorbing "^He on a variety of preplated graphite 
substrates li"— Heat capacity and magnetization measure- 
ments show that the system displays a nearly perfect 
Fermi liquid behavior, with the effective mass to* and 
the enhancement of the spin susceptibility x/xo increas- 
ing with the density. These observations, consistent with 
a divergence of to* near the freezing density, have been 
interpreted- as a signature of a Mott transition leading 
to an insulating crystal. However theoretical approaches^ 
suggest that the singularity of to* and freezing could not 
have the same origin, and indeed the freezing density (as 
well as the magnetic properties of the solid^, and even 
the presence of a possible intermediate phase of a uni- 
form gas of vacancies^) is influenced by the preplated 
substrate. In order to characterize the sole effect of cor- 
relations, it is therefore of particular interest to study 
the ideal, strictly two-dimensional liquid on the verge of 
crystallization, in the absence of any external potential. 
The measured properties of the fluid phase, on the other 
hand, appear to be largely independent of the substrate, 
so that they can be directly compared to the calculated 
properties of the ideal modeli. 

From the theoretical side, such a system provides a 
severe test case for microscopic calculations^, because of 
the strong correlations attained at high densities. We 
thus resort to quantum Monte Carlo (QMC) simulation. 



a powerful tool to study strongly interacting systems, and 
we calculate the ground-state energy per particle & — 
of the 2d ^He liquid at zero temperature as a function of 
the number density p and the spin polarization C,. 

The dependence of the energy on the spin polarization 
is in general very weak in strongly correlated fluids^rJ^. 
The so-called fixed-node (FN) approximationi^, used in 
most QMC calculations, has been argued to give a signif- 
icant bias in the polarization energy of three-dimensional 
liquid ^He^° at high density. 

We thus perform our study going also beyond the 
FN level, following a formally exact methodic,, slightly 
different from the well known transient estimate (TE) 
techniquei^. Briefly, we perform simulations relying on 
the basic Hamiltonian in an enlarged, unphysical space of 
states of any symmetry, including those with Fermi and 
Bose statistics. The ground state energy of the physical 
fermionic '^He is considered as an excitation energy of the 
absolute bosonic ground state, which is sampled exactly 
with QMC. In this approach one trades the sign prob- 
lem faced by TE^"* for the analytic continuation needed 
to extract excitation energies from suitable imaginary- 
time correlation functions. A mixed approach, devised 
to ease detection of the asymptotic convergence of TE by 
a Bayesian analysis of imaginary-time correlation func- 
tions, was proposed By Caffarel and Ceperleyi^. 

In fact a previous FN QMC calculation existsi^, but 
it is limited to low densities and only considers the 
paramagnetic fluid phase. In particular, the accuracy 
of the FN approximation in the high density regime is 
questionable"'^'^. 

We find indeed that the FN level of the theory and 
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the exact calculation predict a qualitatively difTerent be- 
havior: within FN the system becomes ferromagnetic 
well before crystallization takes place upon increasing the 
density, whereas the unbiased calculation shows that the 
spin polarization of the fluid is preempted by freezing, as 
observed experimentally. From the estimated curve e(C) 
we obtain a spin susceptibility enhancement in quantita- 
tive agreement with the available measurements. 



II. QMC SIMULATION 

We simulate A'' particles with the mass ma of ^He 
atoms, interacting with the HFDHE2 pair potential^ in 
periodic boundary conditions. The simulation box, of 
area fi, is a square of side L for the liquid phase; for the 
solid it is a rectangle which accommodates a triangular 
lattice. The Hamiltonian is 
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If the particles obey Bose statistics, projection QMC 
methodsi^"— provide unbiased estimates of the ground- 
state energy and other physical observables. This is made 
possible by the formal similarity between Schrodinger 
equation in imaginary time and the differential equa- 
tion governing a diffusion process in probability theory. 
Such a useful possibility however drops down when Fermi 
statistics enters the game because of the well known sign 
problemi^. 



A. Fixed— Node approach 

The most commonly used approach in the QMC si- 
mulation of Fermions is the fixed-node approximatioiJ^, 
which stochastically solves the imaginary-time 
Schrodinger equation subject to the boundary con- 
ditions implied by the nodal structure of a given trial 
function 'i'T- This approach gives a rigorous upper 
bound to the ground state energy, which often turns out 
to be extremely accurate. 

For a given spin polarization, i.e. considering N-^ spin- 
up and Ni = N — N^ spin-down ^He atoms, ^'t is chosen 
of the form 



= (ri, ...,rjv), T^t = (ri,...,rjv^), 7^_^. 



(2) 



where TZ 

(rTv^+i, ttv), and the whole dependence on the spin 
degrees of freedom is contained in a spin eigenfunc- 
tion for the given polarization 
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describes pair correlations arising from the interaction 
potential; we use a simple pseudopotential of the McMil- 
lan form u{r) = (b/r)^. Finally, the simplest form of the 
antisymmetric factors V {TZf.i) is in the form of Slater 
Determinants of plane waves: 



V (7et4) = det 



^|exp(ifci • rj)| 



(5) 



More accuracy in the FN results is achieved by introduc- 
ing also backflow correlations^^ via quasi-particles vector 
positions: 



V (7et4) = det 



I exp{iki ■ Xj )\ 



(6) 



def 



For the backflow correlation function 77 (r) we adopt the 
simple form: 



(7) 



We will refer to the two choices respectively as plane 
waves fixed-node (PW-FN) and backflow fixed-node 
(BF-FN). For each density, the variational parameters 
6, A, B and C are optimized using correlated sampling22 
at ^ = 0, and left unchanged at different polarizations. 

Part of the bias related to the finite size of the sim- 
ulated system arises from shell effects in the filling of 
single-particle orbitals^'^. This bias can be substan- 
tially reduced adopting twisted boundary conditions^-, 
i.e. choosing k appearing in (O and ([6]) inside the set: 
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where n is an integer vector, L is the side of the simula- 
tion box il and is a twist parameter 9i £ [0, tt] which, 
at the end of the calculations, is averaged over. 

In the solid phase, quantum exchanges are strongly 
suppressed and the energy difference between a Fermionic 
and a Bosonic crystal is negligibly small for the purpose 
of locating the liquid-solid transition. We will there- 
fore replace the energy of ■^He with that of a fictitious 
bosonic Helium of mass ms, which can be calculated 
exactlyi2i22ii^. The small error incurred by such replace- 
ment is bound by the difference between the fermionic 
Fixed-Node (FN) energy and the unbiased bosonic en- 
ergy. This difference, calculated^^ as a check at the melt- 
ing density where it is expected to be largest, is indeed 
in the sub-milliKelvin range. 



B. Fermionic correlations approach 

For the fluid phases the FN approximation may not be 
accurate enough, particularly at high density where cor- 
relations are stronger and the energy balance between dif- 
ferent polarization states is more delicate. Indeed, a FN 
study of three-dimensional '^He, despite the use of sophis- 
ticated backflow correlations, strongly suggests that this 
is the caseiS. We thus perform calculations beyond the 
FN approximation, using a techniqueJ^ which is in prin- 
ciple exact, albeit practically limited to moderate system 
sizes. The idea, in part related to the transient estimate 
formalismpi^ii^, is that of formally viewing ([l} as an oper- 
ator acting inside the Hilbert space 7i{N) = (L^(ri))'*^, 
that is forgetting spin and statistics: one can use Quan- 
tum Monte Carlo to sample the lowest energy eigenfuc- 
tion ^o(TZ) of H among the states of any symmetry. 

It is known^^ that ipo must share the Bose symmetry 
of the Hamiltonian, so that: 



^0 = 



(V'oIV'o)h(A') 



(9) 



is the Ground State energy of a fictitious system of N 
Bosons of mass interacting via the potential v(r). 

The bridge that gives access to fermionic energies may 
be built up as follows: let us fix a spin polarization which 
is surely a good quantum number since the basic Hamil- 
tonian is spin-independent. As discussed in Ref 12, if 
we are able to define an operator Af such that, inside 

nN), 



4'F (7^) = [Apibo) m 



(10) 



has non-zero overlap with the configurational part of any 
exact fermionic Ground State of H for the given ^, we 
can define the imaginary-time correlation function: 



Cf{t) 



,-tH 



Af'>Po)'H{N) 



(V'olV'o)-H(Ar) 



r > (11) 



which can be straightforwardly evaluated in a bosonic 
QMC simulationi^i22ii2i. The lowest energy contribution 
in Cf{t) provides the exact gap between the fermionic 
and the bosonic ground states of the two-dimensional 
Fermi liquid; this can be readily seen by formally ex- 
pressing (fTTj) on the basis {'0n}n>o of eigenvectors of H 
corresponding to the eigenvalues {En}n>Q' 



CF{r) ^ye-(^-<) '<-^f,^°'^"^'' (12) 

A quite natural choiceiS. is to define Af borrowing sug- 
gestions from the form of the trial wave function for the 
FN calculation, i.e.: 
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Figure 1: (Color online) Upper panel: Imaginary time cor- 
relation functions, Cf(t), corresponding to the two different 
choices of the operator in (|13|l . Lower panel: reconstructed 
spectral functions spiijj) obtained with the GIFT method. 



(if Vo) m p(7^t)I?(7^^)Vo(7^) (i3) 

where we can choose either the definition ([SJ of 2? or 
the definition ([6]). We will refer to such choices simply 
as the plane waves fermionic correlations (PW-FC) and 
the backflow fermionic correlations (BF-FC). Naturally 
the flnal results for the Bose-Fermi gap should coincide 
within statistical uncertainties, and the actual compari- 
son can be a good test for the robustness of the approach. 



III. ANALYTIC CONTINUATION 

Once we have achieved a QMC evaluation of Cf (t) , the 
information about the Bose-Fermi gap Abf — Eq—Eq is 
contained in the resulting correlation functions. The re- 
sults for Cf{t) appear as simple smooth decreasing func- 
tions, whose values can be evaluated only in correspon- 
dence with a finite number of imaginary-time values, say 
{tq, ti, T2, t/}; moreover the data are perturbed by un- 
avoidable statistical uncertainties. The Bose-Fermi gap 
Abf thus hidden inside the sets of limited and noisy 
data. How can we extract it? In the upper panel of 
Fig. [1] we show two imaginary time correlation functions 
Cf{t), respectively a PW-FC and a BF-FC, correspond- 
ing to the same spin polarization and twist parameter. 
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The long-T tails of the two curves tend towards a Hn- 
ear behavior (in logarithmic scale) with the same slope, 
and this is a general feature shared by all the functions 
we have evaluated. This indicates that, because of the 
finite-size of the system (and selection rules on the to- 
tal momentum) , the fermionic spectrum has a significant 
gap, i.e. the lowest-energy term exp(— A^f''') in the cor- 
relation function (jl2p appears to be quite well resolved 
with respect to contributions from higher fermionic en- 
ergies. The difference between the two curves (in partic- 
ular the rigid shift between their asymptotic tails) arises 
from the spectral weight of the Ground State contribu- 
tion, which is higher when backflow correlations are taken 
into account, as expected. 

In this favorable situation, the Bose-Fermi gap can be 
reliably extracted by simply fitting an exponential to the 
long-time tail of the correlation function. 

This key result is strongly supported by a more sophis- 
ticated approach, which evaluates IS.bf by performing 
the full Laplace transform inversion of Cf{t), i.e. solving 

p + OO 

Cf{t) = / duje-^'^SFiuj) , (14) 
Jo 

for the unknown spectral function SF{i^)- Recently a new 
method, the genetic inversion via falsification of theo- 
ries (GIFT) method, has been developed to face gen- 
eral inverse problems and in particular it has allowed 
to reconstruct the excitation spectrum of superfluid ^He 
starting from QMC evaluations of the intermediate scat- 
tering function in imaginary-time^; the results were in 
close agreement with experimental data^i. Moreover the 
method has allowed to extract also multiphonon energies 
with a good accuracy level. When applied to the two 
curves depicted in the upper panel of Fig. [1] we find the 
two spectral functions in the lower panel of Fig. [TJ it is 
apparent that the lowest-w peak is indeed well resolved 
from higher-energy contributions. Crucially, its position 
does not depend on the actual choice of the operator 
Af, and it is in excellent agreement with the smallest 
decay constant found by the simple exponential fit. The 
spectral weight instead is different, consistently with the 
differences between PW-FC and BF-FC. 

In this work we adopt an implementation of the in- 
version via falsification of theories, detailed in the Ap- 
pendix, which avoids the rather CPU-intensive genetic 
algorithms'^. This is crucial in the present study, which 
involves an extremely large number of reconstructed 
spectra. Indeed, a single QMC simulation for a given 
density produces correlation functions pertaining to PW 
or BF operators, several spin polarizations, and 15 twist 
parameters in the irreducible wedge of the Brillouin zone 
of the simulation cell; on top of this, data are collected 
in several blocks, individually processed to obtain statis- 
tical uncertainties on the position of the lowest-energy 
peak. 



Table I: The equations of state of '^He for the paramagnetic 
fluid and the solid (solid lines in Figure [2]) are of the form 
aip + a2p^ + asp^ + a^p'^ -I- qs/o^. This Table lists the values 
of the parameters q^, with lengths in A. 





liquid 


sohd 




21.23782 


57.35474 




-1344.413 


-2598.784 


as 


45093.37 


58695.29 




-569306.0 


-532201.7 


as 


4383507 


3063129 
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Figure 2: Equation of state of ^He in two dimensions. Sohd 
line (broken across the coexistence region): liquid and sohd 
''He; dashed line: mass-3 boson fluid; dotted line: liquid '^He, 
after Ref. [l^ . The latter is only reliable at low densities. 



IV. RESULTS 

We fit a fifth order polynomial to the density depen- 
dence of the energies of the triangular crystal and of the 
paramagnetic and the ferromagnetic fluids, listed in Ta- 
ble ID The resulting equation of state of two-dimensional 
■^He is shown in Figured] With the fermionic correlations 
method, we find a transition between the paramagnetic 
fiuid and the triangular crystal around p = 0.061 A~', 
with a narrow coexistence of about 0.002 A~', while the 
ferromagnetic fluid is never stable (see Table |n| . The en- 
ergy of the bosonic mass-3 liquid is also shown. This fic- 
titious system, simulated to extract the PW-FC and BF- 
FC energies, crystallizes at p = 0.069 A~'. The freezing 
density of '^He is considerably higher than the highest 
density simulated in Reffll. Correspondingly, the equa- 
tion of state given in Ref j 16 is only reliable at relatively 
low density. In particular, while it is only slightly below 
our results for p < 0.045 as a consequence of the dif- 
ference of intcrparticle potential adopted'-S, it becomes 
(unphysically) even lower than the bosonic equation of 
state near the melting density, by an amount far larger 
than what could be due to the potential. 

The treatment of the spin polarization state requires a 
special career—. In contrast to Ref.l6l. we find that the 
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Table II: Ground state energy of '^He in K, calculated by the 
FC method for the fluid phases and assumed to equal the 



bosonic energy for the solid phase. 




liquid C = 


liquid C = 1 


solid 


0.020 


0.1707(18) 


0.3218(25) 




0.045 


0.8168(86) 


0.9075(86) 




0.050 


1.1500(81) 


1.2123(88) 




0.055 


1.5972(93) 


1.6574(91) 




0.060 


2.2069(68) 


2.2493(54) 


2.2506(54) 


0.065 


3.0065(73) 


3.0359(45) 


2.9195(26) 


0.070 


4.0644(33) 


4.0915(34) 


3.7878(35) 


0.075 






4.8728(44) 


0.080 






6.2445(35) 


0.085 






7.9589(39) 


0.090 






10.0661(46) 


0.095 






12.6739(39) 


0.100 






15.8536(45) 



BF-FN energy can be significantly higher than the unbi- 
ased Fermionic correlations (FC) energy. Starting from 
negligible values at low density, the BF-FN error quickly 
increases approaching the strongly correlated regime. As 
expected^", it is larger for the paramagnetic than for the 
ferromagnetic fluid. These findings are exemplified in 
Figure [31 The inadequacy of the BF-FN is striking in 
the phase diagram: Figure 2] shows that BF-FN incor- 
rectly predicts a transition to a ferromagnetic fluid well 
before crystallization takes place. Such a transition is 
also evident from Figure [SJ which shows the BF-FN re- 
sults for the polarization energy e{() at various densities. 
The unbiased results, shown in Figure |6l display instead 
a paramagnetic behavior even in a metastable fluid phase 
well beyond the freezing density. 

From the FC polarization energy e(^) we can estimate 
the spin susceptibility enhancement x/xo- Assuming a 
quadratic dependence over the whole polarization range, 
which is generally consistent with the data of Figure [HI 
wc find an excellent agreement with the measured sus- 
ceptibility. Figure [7] shows the comparison between the 
calculated x/Xo s-nd the experimental data. We display 
only the results obtained in the second layer of ^He on 
graphite- since they extend to the highest density in the 
fluid phase, but experiments carried on with differently 
preplated substrates lead to equivalent results in their 
respective density ranges. The agreement among the re- 
sults obtained using different substrates induces us to ex- 
pect that our ideal model actually captures the physical 
mechanisms underlying the behavior of x/xo- 



CONCLUSIONS 
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Figure 3: (Color online) Upper panel: Bose-Fermi gap, 
Abf, as a function of the spin polarization, (, at density 
p =0.020 evaluated via PW-FN, BF-FN, and BF-FC 
with = 18 particles. Middle panel: Bose-Fermi gap, 
Abf, as a function of the spin polarization, (, at density 
p =0.045 A-2 evaluated via PW-FN, BF-FN, and BF-FC 
with A'^ = 18 particles. Lower panel: Bose-Fermi gap, 
Abf, as a function of the spin polarization, at density 
p =0.070 A-2 evaluated via PW-FN, BF-FN, and BF-FC 
with A'' = 26 particles. 

The statistical uncertainties are below the symbols size. 



about 0.002 A^^; the ferromagnetic fluid is never sta- 
ble. From the polarization energy we obtain a spin sus- 
ceptibility enhancement in excellent agreement with the 
experimental values. 

The need for an exact QMC approach is witnessed by 
the failure of the Fixed Node approximation with back- 
flow correlations to predict the lack of a polarization 
transition experimentally observed in the fluid phase, let 
alone an accurate value for the spin susceptibility. 

The estimation of the Bose-Fermi gap via the 
Fermionic correlation method is limited to relatively 
small systems: the present results are obtained with ei- 
ther 18 or (in most cases) 26 particles. While the size ef- 
fect remains the main source of uncertainty of the present 
calculation, the agreement of the calculated and mea- 
sured spin susceptibility suggests that finite-size errors 
are relatively small. 



We have calculated the equation of state and the po- 
larization energy of '^He in two dimensions by means of 
an unbiased QMC method. The system crystallizes into 
a triangular lattice from the paramagnetic fluid at a den- 
sity of 0.061 with a narrow coexistence region of 
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Figure 4: Unbiased FC versus Fixed-Node equation of state. 
Thick solid line (broken across the coexistence region): para- 
magnetic liquid and solid ^He (FC); dashed line: paramag- 
netic liquid (FN); dotted line: ferromagnetic liquid (FN); the 
dashed and dotted lines terminate at the FN freezing density; 
thin solid line: energy of the solid, down to the FN melting 
density. For each density, the energy is relative to the energy 
Eo of the mass-3 boson fluid. 
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Figure 6: (Color online) Exact results for the polarization 
energy -B((") — i?flt(0) relative to the Fermi energy Ep at 
p = 0.020 (open triangles), 0.045 (open squares), 0.050 (filled 
squares), 0.055 (open diamonds), 0.060 (filled diamonds), 
0.065 (open circles), 0.070 (filled circles) A^^, in order of de- 
creasing dispersion. The function i5flt(C) is a quadratic poly- 
nomial in ("^ fitted to the simulation data; the solid line is the 
density-independent result for non-interacting particles. 
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Figure 5: (Color online) Fixed-node results for the polar- 
ization energy E{() — £'flt(0) relative to the Fermi energy 
Ef at p = 0.020 (open triangles), 0.045 (open squares), 
0.050 (filled squares), 0.055 (open diamonds), 0.060 (filled 
diamonds), 0.065 (open circles), 0.070 (filled circles) A~^, i.e. 
from top to bottom. The function Eat (C) is a quadratic poly- 
nomial in fitted to the simulation data; the solid line is the 
density-independent result for non-interacting particles. 



tiative (Laboratory for Interdisciplinary Advanced Simu- 
lation) 2010 grant [http: / /lisa.cilea.it^ 

Appendix A: GIFT algorithm variant 

The inversion procedure that has been employed in 
this work is a variant of the GIFT algorithm^^. This 
new algorithm puts together the idea of the falsification 
principle^^ and a modified implementation of the Prony's 
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Figure 7: Enhancement of the spin susceptibility as a func- 
tion of the density: (filled circles) as measured in the second 
layer of ^He on graphite^; (open circles) as calculated assum- 
ing a quadratic dispersion over the whole polarization range 
in Fig. [6] The corresponding Fixed-node result from Fig. [5] 
would diverge at p ~ 0.050 A~^. 

method^S, a non-iterative parametric technique for mod- 
eling using a linear combination of exponential functions. 
Starting from the basic relation in (fT2|) . which has the 
general form 

oo 

/(T)=^s,e— , (Al) 

1=0 

provided that we are allowed to truncate the previous se- 
ries, X^i^o SfcTo^' Prony's method is computation- 
ally very efficient (it runs in polynomial time) in deducing 
the coefficients {si}"~Q and {u}i}"ZQ from a limited set of 
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estimations, /{kSr) = at fc = 0, 1, . . . , 2n — 1 of /(t), 
being 6t the time step of the QMC simulation. The main 
steps of this algorithm are the following: 



1. solve the regularized linear system 
Ka = b 



(A2) 



defined by the Henkel matrix Kij = f*_^j and by 
the coefficients bi = /*_,_. < n) 

2. find the roots {ziYl^Q ^^'^ polynomial 

z" + a„_iz"-iH haiz + ao (A3) 

as eigenvalues of its respective companion matrix 
and obtain uji — — ^ In 



3. solve the regularized linear system 
As = c 



defined by the Vandermonde matrix Ai 
by the coefficients Ci — f* < n) 



and 



(A4) 



The transition from a nonlinear problem to two linear 
problems and one eigenvalues problem is the main char- 
acteristic of this algorithm; from a mathematical and 
computational point of view this is an advantage. Of 
course, the ill-posedness of this problem remains (our 
implementation uses the truncated singular value decom- 
position regularization) and some care is necessary to 
deal with instability against noised. Such method fits 
the general scheme of the GIFT approach^^, providing a 
very efficient alternative to genetic algorithms in the im- 
plementation of the falsification principle (when dealing 
with Laplace transform inversion). 
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